/****************************************************************************
*
*					Functions of class ExtendedPGMStruct
*
****************************************************************************/

#include <string>
#include <iostream>
#include <fstream>
#include <vector>

#include "ExtendedPGMStruct.h"
#include "SequencingHeuristicFactory.h"
#include "FeatureFactory.h"

#define	minVectorResize		15
// initial value of each weight (is set to be log of the following value)
#define	initWeight			2.0

/****************************************************************************
*					ExtendedPGMStruct::ExtendedPGMStruct
*
*	Description:	Standard constructor with initialization of internal variables
*	Parameters:		-
*	Returns:		-
*
****************************************************************************/
ExtendedPGMStruct::ExtendedPGMStruct (void) {
	sequencingHeuristic = NULL;
}


/****************************************************************************
*					ExtendedPGMStruct::~ExtendedPGMStruct
*
*	Description:	Standard destructor with removing all dynamic internal variables
*	Parameters:		-
*	Returns:		-
*
****************************************************************************/
ExtendedPGMStruct::~ExtendedPGMStruct (void) {
	if (sequencingHeuristic != NULL)
		delete sequencingHeuristic;
}


/****************************************************************************
*					ExtendedPGMStruct::initStructure
*
*	Description:	Initialize internal variables by observing data set parameters
*	Parameters:		environment - a set of parameters from config file
*	Returns:		0
*
****************************************************************************/
int ExtendedPGMStruct::initStructure (Environment& environment) {
	// setting SequencingHeuristic class object
	std::string sqName;
	if (environment.getStringParamValue(sqName, "SequencingHeuristic") != 0) {
		std::cout << "Error in ExtendedPGMStruct::initStructure: cannot find 'SequencingHeuristic' in config file.\n";
		return -1;
	}
	SequencingHeuristicFactory shf;
	sequencingHeuristic = shf.createObjectByName(sqName);
	sequencingHeuristic->initStructure(environment);
	// creating a set of features 
	std::set<std::string> featNames;
	if (environment.getStringSetParamValue(featNames, "Feature") != 0) {
		std::cout << "Error in ExtendedPGMStruct::initStructure: no 'Feature' is set - there is no ingredients to use.\n";
		return -1;
	}
	featureTypes.resize(featNames.size());
	FeatureFactory ff;
	std::set<std::string>::const_iterator it;
	unsigned int i = 0;
	for (it = featNames.begin(); it != featNames.end(); it++) {
		featureTypes[i] = ff.createObjectByName(*it);
		i++;
	}
	// creating a corresponding set of weights
	weights.resize(featureTypes.size(), std::log((float) initWeight));
	return 0;
}


/****************************************************************************
*					ExtendedPGMStruct::initStructure
*
*	Description:	Initialize internal variables by observing data set parameters
*	Parameters:		dataSetMain - data set to initilize internal variables from data
*					dataSetAttr - data set about actors attributes
*	Returns:		0
*
****************************************************************************/
int ExtendedPGMStruct::initStructure (DataSet& dataSetMain, DataSet& dataSetAttr) {
	unsigned int i, j, k;
	// setting data about attributes
	unsigned int numAttr = dataSetAttr.data.size();
	if (numAttr > 0) {
		// using small addition in order to avoid wrong alignment to unsigned int
		nActors = (unsigned int) (std::sqrt ((float) dataSetAttr.data[0].size()) + 0.001);
		actors.resize(nActors);
		for (i = 0; i < nActors; i++) {
			actors[i].attributes.resize(numAttr);
			for (j = 0; j < numAttr; j++) {
				actors[i].attributes[j].resize(nActors, 0);
				for (k = 0; k < nActors; k++) {
					actors[i].attributes[j][k] = dataSetAttr.data[j][i + nActors*k];
				}
			}
		}
	}
	// setting variables and nVars
	setVarValuesFromData (dataSetMain, 0);
	return 0;
}


/****************************************************************************
*					ExtendedPGMStruct::setVarValuesFromData
*
*	Description:	Sets values of variables from a given timestamp of the data set
*	Parameters:		dataSet - data set to initilize internal variables from data
*					timestamp - a timestamp of variables in the data set to use
*	Returns:		0
*
****************************************************************************/
int ExtendedPGMStruct::setVarValuesFromData (DataSet& dataSet, const unsigned int& timestamp) {
	unsigned int i, j;
	nVars = dataSet.maxValues.size() - nActors * (nActors - 1) / 2;
	rates.resize(nVars);
	variables.resize (nVars, 0);
	// setting veriables using data from the first timestamp
	unsigned int nActors2 = nActors*nActors;
	if (timestamp < dataSet.data.size()) {
		// setting Y values
		for (i = 0; i < nActors2; i++)
			variables[i] = dataSet.data[0][i];
		// setting X values
		if (indexOfX.size() == nActors2) {
			// No need to redefine indexOfX
			for (j = 0; j < nActors; j++) {
				for (i = j; i < nActors; i++) {
					variables[getVarID (true, i, j)] = dataSet.data[timestamp][nActors2 + nActors*j + i];
				}
			}
		} else {
			// We need to set indexOfX
			indexOfX.resize( nActors2 );
			for (j = 0; j < nActors; j++) {
				for (i = j; i < nActors; i++) {
					variables[getVarID (true, i, j)] = dataSet.data[timestamp][nActors2 + nActors*j + i];
					indexOfX[getVarID (true, i, j) - nActors2].first = i; 
					indexOfX[getVarID (true, i, j) - nActors2].second = j; 
				}
			}
		}
	} else {
		std::cout << "Error in ExtendedPGMStruct::setVarValuesFromData: there is no timestamp '" << timestamp << "' in the data set.\n";
		return -1;
	}
	// setting number of values of each variable
	numOfValues.resize(nVars);
	for (i = 0; i < nVars; i++)
		numOfValues[i] = dataSet.maxValues[i] + 1;
	// updating info about actors' attributes
	for (i = 0; i < nActors; i++) {
		actors[i].indegree = 0;
		actors[i].outdegree = 0;
		actors[i].degree = 0;
	}
	bool isX; 
	unsigned int iPos, jPos;
	for (i = 0; i < nVars; i++) {
		if (variables[i] == 2) {
			getVarIndex(isX, iPos, jPos, i);
			if (isX) {
				actors[iPos].degree++;
				actors[jPos].degree++;
			} else {
				actors[iPos].indegree++;
				actors[jPos].outdegree++;
			}
		}
	}

	return 0;
}


/****************************************************************************
*					ExtendedPGMStruct::computeRate
*
*	Description:	Returns rate (lambda) for varID and sets all differences 
*					in ingredients for a local net of varID when varID change 
*					its value and when it is does not
*	Parameters:		varID - if of a variable for which to calculate rate
*					ingredients - a vector to fill with differences in 
*						ingredient values
*	Returns:		float - rate value
*
****************************************************************************/
float ExtendedPGMStruct::computeRate (const unsigned int& varID, std::vector <float>& ingredientDiffs) {
	if (varID==1250)
		int j = 0;
	ingredientDiffs.resize(featureTypes.size(), 0.0);
	// if changing variable is not possible at the moment (because of other variables)
	if (!checkVarChange(varID)) {
		for (unsigned int i = 0; i < ingredientDiffs.size(); i++)
			ingredientDiffs[i] = 0.0;
		return 0.0;
	}
	// changing value of the variable
	changeVarValue(varID);
	for (unsigned int i = 0; i < featureTypes.size(); i++)
		ingredientDiffs[i] = featureTypes[i]->calculate (*this, varID);
	// changing the value of the variable back
	changeVarValue(varID);
	for (unsigned int i = 0; i < featureTypes.size(); i++)
		ingredientDiffs[i] -= featureTypes[i]->calculate (*this, varID);
	return computeRate(ingredientDiffs);
}


/****************************************************************************
*					ExtendedPGMStruct::computeRate
*
*	Description:	Returns rate (lambda) for varID (for a local net of varID 
*					when varID change its value and when it is does not)
*	Parameters:		varID - if of a variable for which to calculate rate
*	Returns:		float - rate value
*
****************************************************************************/
float ExtendedPGMStruct::computeRate (const unsigned int& varID) {
	float wsum = 0.0;
	// if changing variable is not possible at the moment (because of other variables)
	if (!checkVarChange(varID)) {
		return 0.0;
	}
	// changing value of a variable
	changeVarValue(varID);
	for (unsigned int i = 0; i < featureTypes.size(); i++)
		wsum += weights[i] * featureTypes[i]->calculate (*this, varID);
	// changing the value of the variable back
	changeVarValue(varID);
	for (unsigned int i = 0; i < featureTypes.size(); i++)
		wsum -= weights[i] * featureTypes[i]->calculate (*this, varID);
	return computeRate(wsum);
}


/****************************************************************************
*					ExtendedPGMStruct::computeRates
*
*	Description:	Initializes rate values for all PGM variables
*	Parameters:		-
*	Returns:		-
*
****************************************************************************/
int ExtendedPGMStruct::computeRates() {
	// Compute all values of vector "rates" defined within the class
	// There are as many rates as there are variables
	for (unsigned int i = 0; i < variables.size(); i++) {
		rates[i].first = computeRate(i, rates[i].second);
		if ( ( rates[i].first > 1.1) && (rates[i].first < 4.9) )
			int jhjh = 0;
	}
	return 1;
}


/****************************************************************************
*					ExtendedPGMStruct::reviseRates
*
*	Description:	Updates rates for all variables after a change in a single variable
*	Parameters:		-
*	Returns:		1
*
****************************************************************************/
int ExtendedPGMStruct::reviseRates (std::map< unsigned int, std::pair<float, std::vector<float>> >& revisedRates, std::vector<unsigned int>& varIDs) {
	unsigned int i, j, k, newVarID;
	bool isX;

	getVarIndex(isX, i, j, varIDs[0]);
	if (i == j) 
		return 1;

	// Update variable value
	changeVarValue (varIDs[0]);
	// Update attributes
	updateAttributes(varIDs[0]);

	// Compute and return rates individually for the variables that are affected by the changes 
	revisedRates.clear();
	
	// Take care of the neighbors of i
	for (k = 0; k < nActors; k++) {
		if ( (k != i) && (k != j) ) {
			// Updating the rates of Y variables
			newVarID = getVarID(false, i, k);
			revisedRates[newVarID].first = computeRate(newVarID, revisedRates[newVarID].second);
			newVarID = getVarID(false, k, i);
			revisedRates[newVarID].first = computeRate(newVarID, revisedRates[newVarID].second);
			// Updating the rates of X variables
			newVarID = getVarID(true, i, k);
			revisedRates[newVarID].first = computeRate(newVarID, revisedRates[newVarID].second);
			//secVarID = getVarID(true, k, i);
			//revisedRates[secVarID] = revisedRates[newVarID];
		}
	}
	// Take care of the neighbors of j
	for (k = 0; k < nActors; k++) {
		if ( (k != i) && (k != j) ) {
			// Updating the rates of Y variables
			newVarID = getVarID(false, j, k);
			revisedRates[newVarID].first = computeRate(newVarID, revisedRates[newVarID].second);
			newVarID = getVarID(false, k, j);
			revisedRates[newVarID].first = computeRate(newVarID, revisedRates[newVarID].second);
			// Updating the rates of X variables
			newVarID = getVarID(true, j, k);
			revisedRates[newVarID].first = computeRate(newVarID, revisedRates[newVarID].second);
			//secVarID = getVarID(true, k, j);
			//revisedRates[secVarID] = revisedRates[newVarID];
		}
	}

	// Take care of the variables on link (i, j)
	// Updating the rates of Y variables
	newVarID = getVarID(false, i, j);
	revisedRates[newVarID].first = computeRate(newVarID, revisedRates[newVarID].second);
	newVarID = getVarID(false, j, i);
	revisedRates[newVarID].first = computeRate(newVarID, revisedRates[newVarID].second);
	// Updating the rates of X variables
	newVarID = getVarID(true, i, j);
	revisedRates[newVarID].first = computeRate(newVarID, revisedRates[newVarID].second);

	// Update variable value
	changeVarValue (varIDs[0]);
	// Update attributes
	updateAttributes(varIDs[0]);
	
	return 1;
}


/****************************************************************************
*					ExtendedPGMStruct::reviseRates
*
*	Description:	Updates rates for all variables after a change in a single variable
*	Parameters:		revisedRates - vector of (variable ID, new rate for variable ID) that 
*						represents variables that change their rate with the change
*						of variable varID
*					revRatesSize - size of the fulfilled 'revisedRates' vector
*					varID - id of a variable that changes its value and forces 
*						other variables to change their rates
*	Returns:		0
*
****************************************************************************/
int ExtendedPGMStruct::reviseRates (std::vector < std::pair <unsigned int, float> >& revisedRates, unsigned int& revRatesSize, const unsigned int& varID) {
	unsigned int i, j, k;
	bool isX;

	getVarIndex(isX, i, j, varID);
	if (i == j) 
		return 1;

	// Update variable value
	changeVarValue (varID);
	// Update attributes
	updateAttributes(varID);

	// Compute and return rates individually for the variables that are affected by the changes 
	revRatesSize = 0;
	
	// Take care of the neighbors of i
	for (k = 0; k < nActors; k++) {
		if ( (k != i) && (k != j) ) {
			// updating size of the revizedRates if needed
			if ((revRatesSize + 3) > revisedRates.size())
				revisedRates.resize(revRatesSize + 3 + minVectorResize);

			// Updating the rates of Y variables
			revisedRates[revRatesSize].first = getVarID(false, i, k);
			revisedRates[revRatesSize].second = computeRate(revisedRates[revRatesSize].first);
			if (rates[revisedRates[revRatesSize].first].first != revisedRates[revRatesSize].second)
				revRatesSize++;
			revisedRates[revRatesSize].first = getVarID(false, k, i);
			revisedRates[revRatesSize].second = computeRate(revisedRates[revRatesSize].first);
			if (rates[revisedRates[revRatesSize].first].first != revisedRates[revRatesSize].second)
				revRatesSize++;
			
			// Updating the rates of X variables
			revisedRates[revRatesSize].first = getVarID(true, i, k);
			revisedRates[revRatesSize].second = computeRate(revisedRates[revRatesSize].first);
			if (rates[revisedRates[revRatesSize].first].first != revisedRates[revRatesSize].second)
				revRatesSize++;
		}
	}
	// Take care of the neighbors of j
	for (k = 0; k < nActors; k++) {
		if ( (k != i) && (k != j) ) {
			// updating size of the revizedRates if needed
			if ((revRatesSize + 3) > revisedRates.size())
				revisedRates.resize(revRatesSize + 3 + minVectorResize);

			// Updating the rates of Y variables
			revisedRates[revRatesSize].first = getVarID(false, j, k);
			revisedRates[revRatesSize].second = computeRate(revisedRates[revRatesSize].first);
			if (rates[revisedRates[revRatesSize].first].first != revisedRates[revRatesSize].second)
				revRatesSize++;
			revisedRates[revRatesSize].first = getVarID(false, k, j);
			revisedRates[revRatesSize].second = computeRate(revisedRates[revRatesSize].first);
			if (rates[revisedRates[revRatesSize].first].first != revisedRates[revRatesSize].second)
				revRatesSize++;

			// Updating the rates of X variables
			revisedRates[revRatesSize].first = getVarID(true, j, k);
			revisedRates[revRatesSize].second = computeRate(revisedRates[revRatesSize].first);
			if (rates[revisedRates[revRatesSize].first].first != revisedRates[revRatesSize].second)
				revRatesSize++;
		}
	}

	if (!isX) {
		// Take care of the variables on link (i, j)
		// updating size of the revizedRates if needed
		if ((revRatesSize + 2) > revisedRates.size())
			revisedRates.resize(revRatesSize + 2 + minVectorResize);
		// Updating the rates of Y variable (do not include Y from i to j - the initial changed ling)
		revisedRates[revRatesSize].first = getVarID(false, j, i);
		revisedRates[revRatesSize].second = computeRate(revisedRates[revRatesSize].first);
		if (rates[revisedRates[revRatesSize].first].first != revisedRates[revRatesSize].second)
			revRatesSize++;
		// Updating the rates of X variables
		revisedRates[revRatesSize].first = getVarID(true, i, j);
		revisedRates[revRatesSize].second = computeRate(revisedRates[revRatesSize].first);
		if (rates[revisedRates[revRatesSize].first].first != revisedRates[revRatesSize].second)
			revRatesSize++;
	}

	// Update variable value
	changeVarValue (varID);
	// Update attributes
	updateAttributes(varID);
	
	return 0;
}


/****************************************************************************
*					ExtendedPGMStruct::updateRates
*
*	Description:	Updates revised rates for all variables enduring changes 
*					after a change in a single variable and puts changes to rates, 
*					stored in this class
*	Parameters:		varIDs - indexes of variables that changes their values
*					TERM2 = sum_{j} lambda_j * (fTheta(PGMe with j) - fTheta(PGMe without j)
*	Returns:		float - difference in sum of rates after and before updates
*
****************************************************************************/
float ExtendedPGMStruct::updateRates (const unsigned int& varID, std::vector <float>& TERM2) {
	unsigned int i, j, k, m, newVarID;
	bool isX;
	float toRet = 0.0;


	// Compute and return rates individually for the variables that are affected by the changes 
	getVarIndex(isX, i, j, varID);
	if (i == j) 
		return 1;

	// Update variable value
	changeVarValue (varID);
	// Update attributes
	updateAttributes(varID);
	
	// Take care of the neighbors of i and j
	for (k = 0; k < nActors; k++) {
		if ( (k != i) && (k != j) ) {
			// Updating the rates of Y variables
			newVarID = getVarID(false, i, k);
			toRet -= rates[newVarID].first;
			for (m = 0; m < weights.size(); m++)
				TERM2[m] -= rates[newVarID].first * rates[newVarID].second[m];
			rates[newVarID].first = computeRate (newVarID, rates[newVarID].second);
			toRet += rates[newVarID].first;
			for (m = 0; m < weights.size(); m++)
				TERM2[m] += rates[newVarID].first * rates[newVarID].second[m];

			newVarID = getVarID(false, j, k);
			toRet -= rates[newVarID].first;
			for (m = 0; m < weights.size(); m++)
				TERM2[m] -= rates[newVarID].first * rates[newVarID].second[m];
			rates[newVarID].first = computeRate (newVarID, rates[newVarID].second);
			toRet += rates[newVarID].first;
			for (m = 0; m < weights.size(); m++)
				TERM2[m] += rates[newVarID].first * rates[newVarID].second[m];

			newVarID = getVarID(false, k, i);
			toRet -= rates[newVarID].first;
			for (m = 0; m < weights.size(); m++)
				TERM2[m] -= rates[newVarID].first * rates[newVarID].second[m];
			rates[newVarID].first = computeRate (newVarID, rates[newVarID].second);
			toRet += rates[newVarID].first;
			for (m = 0; m < weights.size(); m++)
				TERM2[m] += rates[newVarID].first * rates[newVarID].second[m];

			newVarID = getVarID(false, k, j);
			toRet -= rates[newVarID].first;
			for (m = 0; m < weights.size(); m++)
				TERM2[m] -= rates[newVarID].first * rates[newVarID].second[m];
			rates[newVarID].first = computeRate (newVarID, rates[newVarID].second);
			toRet += rates[newVarID].first;
			for (m = 0; m < weights.size(); m++)
				TERM2[m] += rates[newVarID].first * rates[newVarID].second[m];

			// Updating the rates of X variables
			newVarID = getVarID(true, i, k);
			toRet -= rates[newVarID].first;
			for (m = 0; m < weights.size(); m++)
				TERM2[m] -= rates[newVarID].first * rates[newVarID].second[m];
			rates[newVarID].first = computeRate (newVarID, rates[newVarID].second);
			toRet += rates[newVarID].first;
			for (m = 0; m < weights.size(); m++)
				TERM2[m] += rates[newVarID].first * rates[newVarID].second[m];
			
			newVarID = getVarID(true, j, k);
			toRet -= rates[newVarID].first;
			for (m = 0; m < weights.size(); m++)
				TERM2[m] -= rates[newVarID].first * rates[newVarID].second[m];
			rates[newVarID].first = computeRate (newVarID, rates[newVarID].second);
			toRet += rates[newVarID].first;
			for (m = 0; m < weights.size(); m++)
				TERM2[m] += rates[newVarID].first * rates[newVarID].second[m];
		}
	}

	// Take care of the variables on link (i, j)
	// Updating the rates of Y variables
	newVarID = getVarID(false, i, j);
	toRet -= rates[newVarID].first;
	for (m = 0; m < weights.size(); m++)
		TERM2[m] -= rates[newVarID].first * rates[newVarID].second[m];
	rates[newVarID].first = computeRate (newVarID, rates[newVarID].second);
	toRet += rates[newVarID].first;
	for (m = 0; m < weights.size(); m++)
		TERM2[m] += rates[newVarID].first * rates[newVarID].second[m];

	newVarID = getVarID(false, j, i);
	toRet -= rates[newVarID].first;
	for (m = 0; m < weights.size(); m++)
		TERM2[m] -= rates[newVarID].first * rates[newVarID].second[m];
	rates[newVarID].first = computeRate (newVarID, rates[newVarID].second);
	toRet += rates[newVarID].first;
	for (m = 0; m < weights.size(); m++)
		TERM2[m] += rates[newVarID].first * rates[newVarID].second[m];

	// Updating the rates of X variables
	newVarID = getVarID(true, i, j);
	toRet -= rates[newVarID].first;
	for (m = 0; m < weights.size(); m++)
		TERM2[m] -= rates[newVarID].first * rates[newVarID].second[m];
	rates[newVarID].first = computeRate (newVarID, rates[newVarID].second);
	toRet += rates[newVarID].first;
	for (m = 0; m < weights.size(); m++)
		TERM2[m] += rates[newVarID].first * rates[newVarID].second[m];

	return toRet;
}


/****************************************************************************
*					ExtendedPGMStruct::updateAttributes
*
*	Description:	Updates network structure-dependent actor attributes 
*					(indegree, outdegree, degree) after a change in a single variable
*	Parameters:		-
*	Returns:		0
*
****************************************************************************/
int ExtendedPGMStruct::updateAttributes (const unsigned int& varID) {
	unsigned int i, j = 0;
	bool isX;

	getVarIndex(isX, i, j, varID);

	// Update the count
	if (isX) {
		if (variables[varID] == 2) {
			actors[i].degree++;
			actors[j].degree++;
		} else {
			actors[i].degree--;
			actors[j].degree--;
		}
	} else { // VarID - is Yij 
		if (variables[varID] == 2) {
			actors[i].indegree++;
			actors[j].outdegree++;
		} else {
			actors[i].indegree--;
			actors[j].outdegree--;
		}
	}
	return 0;
}


/****************************************************************************
*					ExtendedPGMStruct::clear
*
*	Description:	Removes all variables in object internal variables
*	Parameters:		-
*	Returns:		-
*
****************************************************************************/
void ExtendedPGMStruct::clear (void) {
	PGMStruct::clear();
	nActors = 0;
	actors.clear();
	rates.clear();
	if (sequencingHeuristic != NULL)
		delete sequencingHeuristic;
}


/****************************************************************************
*					ExtendedPGMStruct::learn
*
*	Description:	Main Learning Procedure
*	Parameters:		dataSet - a set of data samples
*	Returns:		0 - if there were no errors during the optimization
*
****************************************************************************/
int ExtendedPGMStruct::learn (DataSet& dataSet) {
	unsigned int j, n;
	int ret = 0;

	// Getting Number of Weights
	n = weights.size();

	// Representing weights in form of lbfgsFLOATval_t*
    lbfgsfloatval_t *x = lbfgs_malloc(n);
	// Checking that we've got the memory
    if (x == NULL) {
        printf("ERROR in ExtendedPGMStruct::learn: Failed to allocate a memory block for weight variables.\n");
        return -1;
    }
    // Initialize the variables
	for (j = 0; j < weights.size(); j++)
		x[j] = (lbfgsfloatval_t) weights[j];

    // Initialize the parameters for the L-BFGS optimization
	lbfgs_parameter_t lbfgsParams;
    lbfgs_parameter_init(&lbfgsParams);
	lbfgsParams.linesearch = 2;
	lbfgsParams.epsilon = 0.01;
	lbfgsParams.ftol = 0.01;
	lbfgsParams.max_linesearch = 5;
	//lbfgsParams.linesearch = LBFGS_LINESEARCH_BACKTRACKING;

	// Returned optimal functional value
	lbfgsfloatval_t fx; 
	RequiredStructures requiredStructures;
	requiredStructures.dataSet = &dataSet;
	requiredStructures.extPGMStruct = this;

	// Start the L-BFGS optimization; this will invoke the callback functions calcValueAndGradient() and printoutProgress() when necessary.
	// This function finds minimum values of fx as a function from x
    ret = lbfgs(n, x, &fx, calcValueAndGradient, printoutProgress, &requiredStructures, &lbfgsParams);

    // Report the result
    /*printf("\n");
	printf("L-BFGS optimization terminated with status code = %d\n", ret);
    printf("  fx = %f\n", fx);
	for (j = 0; j < n; j++)
		printf("  x[%d] = %f", j, x[j]);
	printf("\n");*/

	// Getting parameters back to the matrix 'weights'
	for (j = 0; j < weights.size(); j++) {
		weights[j] = (float) x[j];
	}

	// Freeing the memory of x
    lbfgs_free(x);

	return ret;
}


/****************************************************************************
*					ExtendedPGMStruct::calcValueAndGradient
*
*	Description:	Function calculates values of objective function and its 
*					gradient at a given set of weights values
*	Parameters:		instance - data sent for lbfgs() function (all sturctures)
*					x - current values of variables
*					g - current gradient values of variables
*					n - number of variables
*					step - line-search step used for this iteration
*	Returns:		lbfgsfloatval_t - values of the objective function in this
*						set of parameters (weights)
*
****************************************************************************/
lbfgsfloatval_t ExtendedPGMStruct::calcValueAndGradient (
    void *instance, // our structures to calculate values
    const lbfgsfloatval_t *x, // current vector of feature weights
    lbfgsfloatval_t *g, // gradiend (needs to be filled)
    const int n, // number of parameters
    const lbfgsfloatval_t step ) // current step of the search routine
{
	// casting 'instance' to be of type PGMandDataSet*, as we created it in 'learn' function
	RequiredStructures* requiredStructures;
	requiredStructures = (RequiredStructures*) instance;

	// setting new theta values
	for (int j = 0; j < n; j++)
		requiredStructures->extPGMStruct->weights[j] = (float) x[j];
	
	// calling sequencing function to get probability and gradients
	std::vector<unsigned int> sequence; 
	float seqLogLikelihood; 
	std::vector<float> gradients;
	requiredStructures->extPGMStruct->sequencingHeuristic->findSequence (sequence, seqLogLikelihood, gradients, *(requiredStructures->dataSet), *(requiredStructures->extPGMStruct));
	
	// setting gradient values
    for (int i = 0; i < n; i ++) {
        g[i] = - gradients[i];
    }
    return (lbfgsfloatval_t) - seqLogLikelihood;
}


/****************************************************************************
*					ExtendedPGMStruct::printoutProgress
*
*	Description:	Function printout current state of optimization
*					Specifically it prints current values of the objective 
*					function and weights and convergence parameters
*	Parameters:		instance - data sent for lbfgs() function (all sturctures)
*					x - current values of variables
*					g - current gradient values of variables
*					fx  - current value of the objective function
*					xnorm - the Euclidean norm of the variables
*					gnorm - the Euclidean norm of the gradients
*					step - line-search step used for this iteration
*					n - number of variables
*					k - iteration count
*					ls - number of evaluations called for this iteration
*	Returns:		0
*
****************************************************************************/
int ExtendedPGMStruct::printoutProgress (
    void *instance,
    const lbfgsfloatval_t *x,
    const lbfgsfloatval_t *g,
    const lbfgsfloatval_t fx,
    const lbfgsfloatval_t xnorm,
    const lbfgsfloatval_t gnorm,
    const lbfgsfloatval_t step,
    int n,
    int k,
    int ls )
{
	// Report the result
    printf("\n");
	printf("Iteration %d of optimization:\n", k);
    printf("  fx = %f\n", fx);
	for (int i = 0; i < n; i++)
		printf("  x[%d] = %f,   g[%d] = %f", i, x[i], i, g[i]);
	printf("\n");
	printf("  xnorm = %f, gnorm = %f, step = %f\n", xnorm, gnorm, step);
    return 0;
}
